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Methods and results for numerical simulations of one and two interacting rf-SQUID systems 
suitable for adiabatic quantum gates are presented. These are based on high accuracy numerical 
solutions to the static and time dependent Schroedinger equation for the full SQUID Hamiltonian in 
one and two variables. Among the points examined in the static analysis is the range of validity of 
the effective two-state or "spin 1/2" picture. A range of parameters is determined where the picture 
holds to good accuracy as the energy levels undergo gate manipulations. Some general points are 
presented concerning the relations between device parameters and "good" quantum mechanical state 
spaces. 

The time dependent simulations allow the examination of suitable conditions for adiabatic behav- 
ior, and permits the introduction of a random noise to simulate the effects of decoherence. A formula 
is derived and tested relating the random noise to the decoherence rate. Sensitivity to device and 
operating parameters for the logical gates NOT and CNOT are examined, with particular attention 
to values of the tunnel parameter (3 slightly above one. It appears that with values of [3 close to 
one, a quantum CNOT gate is possible even with rather short decoherence times. 

Many of the methods and results will apply to coupled double-potential well systems in general. 

INTRODUCTION 

In previous work we have described quantum logic gates based on the rf SQUID. The basic operation involved is an 
adiabatic inversion, where the SQUID reverses flux states under the sweep of an external field 0'^^* . This is equivalent 
to the logical NOT ^. When a second SQUID is added whose flux can add or subtract from 0'^^*, parameter ranges 
were found for the two so that the two-SQUID system undergoes a reshuffling of levels equivalent to the 

logical operation CNOT [|,[|. 

In this paper we present a study of these systems by numerical simulations which enable us to examine these 
processes in more quantitative detail. Among the points we can examine is the validity of the "spin 1/2 picture". In 
the previous work we often found good agreement with a simplified "spin 1/2 analogy" where the two lowest states of 
each SQUID are treated as an effective two-state system. This picture is very useful in understanding and predicting 
the behavior of the systems and here we examine its validity by simulations for the full many-state system. 

A further question which can be studied in detail is adiabaticity. This is the operating principle for our quantum 
gates and it is necessary to know under which conditions it holds. 

Finally, we can use our programs to study the effects of decoherence on the quantum gates. We will examine how 
to introduce decoherence as a noise signal and study its effects on our operations. 

It should be stressed that our results would have been difficult if not impossible to obtain without the extensive 
system of numerical programs. Due to the sensitivity to the various parameters and the subtleties of the tunneling 
problem, analytic methods would be difficult and uncertain. With the high accuracy programs, a short run can 
replaces otherwise long, complicated, and often approximate formulas. 
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ONE VARIABLE SYSTEM 

We begin by shortly reviewing our approach (H, 01 applied to the one variable, one rf SQUID system, at first 
without decoherence. This will allow us to fix the notation and parameters, and to establish the connection between 
the full Hamiltonian and the "spin 1/2 picture". 



Squid Hamiltonian 

The SQUID Hamiltonian in terms of the capacitance C and inductance L of the junction is 7i = — ~t5T7^2 + U 

2C( — )- f 

with U = ~ 0"^*)^] + l3cos(j)}, where /3 = 2^^, Ic being the critical current for the junction 4]. The 

variable cj) is the flux <I> in the SQUID loop in flux quantum units (f> — and furthermore shifted by tt so that 

(j) = corresponds to the maximum at the center of the double potential well U . Similarly cf)'^''^^ is an external field, 
where due to the shift, (j)'^^'^ — corresponds to a non-zero applied field of size By factoring out an overall energy 
scale Eq = 1/^/LC so that H = EqH wc obtain the dimensionless hamiltonian 

H = + Vomcb - cj^^^r] +(3coscb} = ^-^+ - </'^^*)'] + /3 . (1) 

With this choice of the factor Eq the "mass" /i and the "potential" Vq are equal, hence the second form in terms of a 
common parameter ^, whose value is (, — fi — Vq ^ 1030 y^^^^. In effect the capacitance C and inductance L have 
been exchanged for an energy scale Eq and a dimensionless number ^. We shall discuss the physical meaning of ^ 
below. 

In the following we shall endeavor to express all energies in terms of the general energy unit Eq and all times in the 
time unit 1/Eq. To converted these to dimensional units one multiplies by Eq = 6.4 x 10~'*ey x [L/pH C IpF)^^!"^ — 
1.0 X 1012 radiansjsec x (L/pff C/pF)-^/^ = 7.7K x (L/pH C/pF)-^/^ and for the time by 1/Eq = 1.0 x 
10~^'^ seconds x [L/pH C/pFY^'^ Thus results with the dimensionless Hamiltonian, Eq[l] which we will use in 
our computer simulations, are to be converted to physical energies by multiplying by Eq. Typical values L — AQQpH 
and C = Q.lpF for example, yield Eq — 160 x XO'^radians/ sec — 1.2 K and ^ = 16.3, while the time unit is 
6.3 X 10~i2 sec. We will use these values for "typical examples" or when absolute times or energies are needed. We 
shall usually work with (3 ~ 1.19, which for L — AOOpH corresponds to a critical current of about 1 iiA. 

To carry out calculations with Eq[l]for the eigenvalues and eigenfunctions, the numerical method employed a large 
basis of harmonic oscillator wavefunctions and expanded the cosine as a low order polynomial. Since in the harmonic 
oscillator basis polynomials are simple, sparse, matrices, the problem is reduced to algebraic manipulations and a 
matrix diagonalization. It was usually found that an expansion up to 8th order for the cosine and a basis of 256 
oscillator states gave numerical stability. A typical run with a few hundred basis vectors on a Pentium 4 machine 
lasts less than a minute and resolves our smallest energy splitting to better than two significant figures. 



Parameters of the Hamiltonian 



Eq[l]contains three parameters, 0^^*, /3 and The applied flux controls the asymmetry of the potential and for 
0ea;t_Q ^]^^ potential is symmetric. Fig. 1, Left, shows an example with the potential in this symmetric configuration. 
In this "level crossing" situation the energy splitting of the lowest pair of levels is just the tunneling energy ujtunnei 
and so is quite small. Here it is 0.0044 and not clearly resolved on the plot. 

Fig. 1, Right, shows, under the same conditions, an asymmetric configuration with (^'^^*=0.0020. To produce an 
adiabatic inversion or NOT, (p^^'^ begins at such a non-zero value and adiabatically "sweeps" to the opposite value. 
The level splitting is now dominated by the shift of the potential and not the tunneling. The splitting of the lowest 
pair is 0.060, an order of magnitude greater than in the symmetric configuration. The wavefunctions of the energy 
eigenstates are now concentrated in the left or right well, while at "level crossing" they were ± linear combinations 
of these states. 

Turning to the parameter ^, raising its value both increases the "mass" and the height of the potential and so will 
generally lead to a greater concentration of the wavefunction in one of the potential wells, as well as a reduction of 
the tunneling. In the wells one has roughly harmonic oscillator potentials, for which the ground state wavefunction is 
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FIG. 1: Potential with first four energy levels according to Eq[T] for {3 — 1.19 and Vb = /i = 16.3. Left: The symmetric or 
"level crossing" configuration with (fi'^^^—Q. The lowest line represents a pair of levels which are not clearly resolved on the plot. 
Right: An asymmetric configuration with (/!>^'^'=0.0020, showing a much greater pair separation. The horizontal axis represents 
the reduced flux the vertical axis the energy in units of 6.4 X lQ-'^eV{L/pH C/pF)-^'^. 



ip e 2'\/^^^<A ^ where 0' means (p referred to its equilibrium position and c is a factor of order one. Thus ^ gives 
the spread of the wavefunction in the well; it controls how "uncertain" the flux is when it is in a "definite" state. 
Similarly the exponent in the tunnehng factor e~ -f '^^ V^/Icvb-^ increases with ^ and reduces the tunneling. Therefore 
^ = \/ fiVo may be viewed as a measure of how "classical" the wavefunction is, and the relatively large value of ^ we 
deal with indicates that our component wavefunctions are rather well defined and "classical" in this sense. That is 
although we deal with linear combinations of states, these basis states themselves are relatively localized, are rather 
"classical" . Examples of this will be presented in Fig. 9. 

The parameter /?, finally, has a strong effect on the tunneling since increasing its value both widens and heightens 
the barrier. Thus increasing f3 from 1.19 to 1.35 leads to three pairs of well defined states below the barrier, with a 
tunnel splitting of only 2.0 x 10^^ for the lowest pair. In the other direction, for /3 somewhat less than one, there are 
no longer any states localized around a definite value of the flux at all. The most interesting region for the present 
purposes appears to be the values of /3 somewhat greater than one, and for most of our simulations we shall use 
(3 = 1.19. With f3 = 1.19, ^ = 16.3, the splitting at level crossing is 0.0044, while the distance to the next set of levels 
is 0.48. This difference of two orders of magnitude provides a factor of ten margin in manipulating the levels with 
still a factor of ten to the next set. However in studying specific designs it may be of interest to make a careful study 
of the behavior with respect to (3 in the vicinity of one, and we shall also briefly consider (3 = 1.14. 

Finally, it should be observed that the energy splitting from our lowest pair of states to the next ones is on the 
order of the energy unit, which is Eq 1 K . Consequently when working well below IK, which we shall assume, one 
can expect only the lowest pair of states to be populated, and the neglect of thermal over-the-barrier transitions to 
be justified. 



IDENTIFICATION WITH A SPIN 1/2 SYSTEM 

It is frequently a useful simplification to use the "spin picture" where we treat two closely separated levels, such as 
the lowest pair in Fig. 1, as the two states of a "spin 1/2 system" and for the two qubit system as two such spins. 
In this section we discuss the relation between this picture and the full Hamiltonian. The spin picture Hamiltonian 
in the absence of noise or decoherence is 

H=^a-Y, (2) 

where we drop an irrelevant additive constant. 

We wish to use the two lowest states of Eq[T]as our qubit, and to identify it with a system which can be effectively 
described by the simple Eq[21 In other words, we wish to use the two lowest energy eigenstates found from the exact 
Eq[l] with say 0'^^*=O, to span a two-dimensional basis to set up the spin picture. We assume that as 0^^^* is varied 
over a small range we stay in this Hilbert space and always have to do with various linear combinations of the same 
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wavefunctions. This seems a plausible assumption when the tunneling energy is small compared to the other level 
splittings, and we shall present evidence for it below. 

Having made this assumption, the next question is which linear combinations are to be used as the fixed basis-in 
the spin language which linear combinations to chose as "spin up" and "spin down" along an abstract "z axis" . As 
for the flavor with neutrinos and K mesons, or the handedness with chiral molecules Q, we wish to chose these basis 
states to be eigenstates of a definite, externally measurable, property of the system. Here we choose the direction of 
the flux, i.e. the current in the SQUID, corresponding to the system localized either in the left or the right potential 
well. Due to the tunneling these states are not energy eigenstates and so not stationary; they will generally undergo 
oscillations in time. 

This choice implies that the basis states, those to be identified with "spin up" and "spin down" along the "z axis" , 
are chosen to be eigenstates of the flux cj). In terms of the full Eq [U these are wavefunctions concentrated in one 
potential well only. Now in our two-dimensional space spanned by the two lowest eigenstates, the "position" variable 
(/) is a 2 X 2 matrix. This matrix can be expanded in terms of the Pauli matrices. Our choice that the basis states are 
eigenstates of means that the abstract axes are defined such that (j) is proportional to the diagonal Pauli matrix tj^ 



~ o-z • (3) 

Furthermore, we can approximately establish the constant in this relation by using the fact that for the so-defined 
eigenstates of cr^ the wavefunction is approximately localized. We then approximately know the value of 0, since 
operating on one of these states, for example \R > for a state concentrated in the right well, we expect >^ (j^clR > 
where the number (j)c is the value of the variable </> where the wavefunction is concentrated, say has its maximum 
value. Since az operating on one of the localized states has been defined as ±1, Eq[3]becomes 

(j) w (f>ccrz (4) 

As may be seen in Fig. 1, 4>c will usually be roughly one. 

With this identification we now proceed to analyze the full Hamiltonian Eq[T] Since we always take 0^^^* small, it 
is a good approximation to write Eq[l]as 

^^+e{|0'+/3W}-C0^^V- (5) 

Eq [5] represents the total Hamiltonian as the Hamiltonian for the case of the symmetric or "level crossing" potential 
where (lf^^=Q and an asymmetric term proportional to (j). Evidently, in view of Eq|3]and the symmetry of the first 
part of Eqini the last term in Eq[5]is to be identified with the ^VzCTz of the spin Hamiltonian Eq[51 

From this observation we can, by using Eq 31 identify the value of Vz in Eq [2] with the parameters of the full 
Hamiltonian Eq[l]as 



Vz - 2e0^"Vc ■ (6) 

With (j)c close to one, we anticipate Vz ~ 2^(/)^^*. Below we show a more exact calculation of Vz from the numerical 
results with the full Hamiltonian. 

The quantity from the full Hamiltonian to be associated with the magnitude V of V in the spin picture is easy to 
identify. The eigenvalue splitting of ^cr • V is always V. Thus given a numerical calculation with Eq [T] yielding the 
level splittings, we anticipate 

splitting = V = \/vf+Vl , (7) 

where splitting is the energy difference of the two lowest states. (A Vy does not enter into our considerations since 
all quantities are real and ay would involve imaginary quantities.) 



The x-component Vx represents the tunneling energy oJtunneU or the splitting resulting from the first, symmetric 
part of Eq [Sj - the Hamiltonian at "level crossing" . We can therefore obtain Vx from the numerical evaluation at 

^'^^'cx Vz = 0: 

Vx = i^tunnel = Splitting\(^,^t^Q) . (8) 

The numerical value of Vz then follows from those for V and Vx ■ 
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TABLE I: Quantities from the exact numerical simulation for comparison with those of the "spin 1/2 picture". The energy 
splittings are to be seen in relation to the distance to the next set of levels, which is about 0.6 units. The angle dv refers to the 
angle the V vector makes with the z-axis, while the angle 0^ refers to the same angle inferred from the spinor eigenfunctions, 
as described in the text. "Completeness" characterizes how well the two states in question remain in the same state space as 
0'^^*is varied. In general, the results indicate that the two state picture holds to about <ff^^Ki 0.01. Blank spaces in tables imply 
repetition of previous values. 

Energy level behavior 

We can test this picture, where Vx and Vz are treated as approximately independent quantities, by plotting the 
splittings from the full numerical calculation vs (f)'^^^ to check if they have the expected form ^/V^ + V^, with Vz 
proportional to (jf^^ and Vx a constant. This is shown in Fig. 2. There is an excellent fit to this form with the fit 
constant for Vx equal to the value at level crossing. In addition the fit coefficient C in Vz =2C(/)'^^* is 14.9, while from 
the estimate EqlHlwe would expect C ~ Vb^c which with (/)c ~ 1 is 16.3. Or if we adjust to make the identification 
exact, we need (j)c = 0.91. Inspection of a plot of the wavefunction shows that its maximum is indeed close to this, at 
about (f) = 0.88. 

Table 1 shows some of these values and also those for some larger cjf^*. Deviations from the linear behavior 
splitting (x(jf^* begin to set in at about (f)^^*= 0.015, where the splitting is near a half unit. As would be expected, 
this is on the order of the distance to the next set of levels, as seen in Fig. 1. Thus, for the behavior of the energy 
levels at small 0*^^* , there is quite good agreement between the numerical calculations with the full Hamiltonian and 
the "spin 1/2 picture". Experimentally, a plot equivalent to Fig. 2 has been mapped out to about 0'^^*=O.OO8 

Rotation angle 

As another check on the simple two-state picture, we can consider two independent ways of finding the angle the 
"spin" makes with the abstract z axis. One way is to use Vx and the magnitude V from the energy splittings. These 
are given by Eqs [8] and [7| and the resulting 9v from sinOy = Vx/V is given in Table I. A second way, using only the 
Hamiltonian at a given 0'^^*, is to find the angle of rotation from the energy eigenstates to the "ctz" eigenstates. If 
the angle that V makes with the z axis is 6, then the eigenstates v± of the spin Hamiltonian are given in terms of the 
spin-up, spin -down states u± as (cos|m+ -|-sm|u_) and (— sin|w+ +cos^U-). Thus by finding the rotation from the 
V to the u we may determine 9. Numerically, this is done by computing the 2x2 matrix of </) in the energy eigenstate 
basis and finding the rotation necessary to diagonalize it. The results are shown in the Table as 6^. Disagreement 
between the two methods begins to appear around (/)'^^*=0.01, where the splitting is 0.3 energy units. This comparison 
is perhaps more sensitive than that using the energy levels alone. However the disagreement is only significant when 
the angle is small. 

Hilbert space completeness 

Finally, as another check, we can try to directly see if the system remains in the same two dimensional Hilbert space 
as c/)'^^* is varied. We examine this by comparing the spaces spanned by the two lowest eigenstates of the Hamiltonians 
with different 0^^^*. The worst overlap between the eigenstates of the Hamiltonian with 0'^^*=O and those with the 
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FIG. 2: Level splitting as a function of (^"^^^ for small 0*^^*, from the numerical calculation with the full Hamiltonian. A good 
fit with the "spin 1/2" form splitting = ^^(^1 + (K)^ = {2C + (14)2 is obtained, yielding K = 0.0040 and C = 14.9. 
C is in agreement with the prediction Eq|6l 

non-zero ^'^^'are listed in the column "Completeness". The method will be explained below when the two dimensional 
case is discussed. Significant deviations again appear around 0'^^*=O.O15 and one notes a sudden change as mixing 
with the next set of principal states becomes important. 

In summary, the two state picture seems to work well with the parameters used here, up to about 0.01 for 0*^^*, or 
pair splitting ~ 0.2 energy units, to be compared with the principal level spitting of ^ 0.6 units. Inside this range, 
the picture that we always deal with different linear combinations of the same two states while the Hamiltonian varies 
seems to be justified. 

ADIABATICITY 

Our gates operate by a sweep of the externally applied (^f^^. The speed of a sweep is of course relevant to how fast 
a device or a set of devices might operate. Probably more important than the simple speed, however, is its connection 
with the decoherence question. The decoherence is characterized by a rate (our D below). Therefore fast gates, giving 
the decoherence less time to act, are favorable from the point of view of decoherence. 

Although in principle very fast sweeps thus seem desirable, this is not possible without violating the adiabatic 
condition upon which the gate operations are based. It is therefore of interest to find out how fast sweeps can be 
performed without violating adiabaticity. Using the simulation programs we can study this point in detail. First we 
examine the simple case of the adiabatic inversion or NOT, using one SQUID without decoherence. 

To deal with the time dependent Hamiltonians numerically, repeated iterations of (1 + iHAt) were employed, 
applied to wavefunctions found by the methods of the static calculation described above. To save time in such runs 
the matrix H was calculated not in the original large oscillator basis but in a reduced basis of the few lowest energy 
eigenstates. The results could be checked by incorporating more states into this "second cut" basis. Usually, as might 
be expected from the arguments around Table I, two states were sufficient in one variable and four states in the two 
variable problem. 

According to the estimate given in ref 0] , adiabaticity is guaranteed when the sweep time t sweep is sufficiently long 
such that 

<< ^tunnel J (9) 

f^sweep 

where e is the initial energy level splitting (same as the initial difference of the minima of the potential wells to 
about 10%). The condition is sensitive to oJtunnei because the violation of adiabaticity takes place essentially at level 
crossing. Violation occurs when the rate of change is too fast; and the left-hand-side of Eq[9] characterizes this velocity 
of change of the system. Note we cannot make the left-hand-side smaller by simply reducing e. If we want to retain a 
good definition of the original state as approximately a state of definite fiux, one requires Vz/Vx >> 1 or e >> cutunnei, 
that is a substantial initial splitting compared to the splitting at "crossing" . 
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FIG. 3: Study of adiabaticity. The behavior of the scalar product (squared) between an evolved wavefunction and the target 
eigenstate of the final Hamiltonian. The evolved wavefunction is calculated from the Schroedinger equation with the Hamiltonian 
Eq[T] containing a time- varying (j)"^*^ that reverses sign. Two sweeps are shown, one (lower) with tg^eep ~ 7500 time units and 
one (upper) with tsweep = 15000 units. "Level crossing" or (^'^^*=0 occurs at ^tsmeep- The adiabatic inversion or reversal of 
the flux state of the SQUID is more "successful" for the slower sweep. Conditions as for the (3 = 1.19 case in Table II. 

We may characterize the "success" of an adiabatic inversion by starting with a wavefunction which is an eigenstate 
of the initial potential, evolving it in time with the changing potential, and finally comparing it with the stationary 
wavefunction to which it should arrive, namely the eigensolution for the final potential. In the case of perfect 
adiabaticity the two wavefunctions should be the same. Thus we can gauge the loss of adiabaticity by the deviation 
from one of the scalar product of the evolved wavefunction with the "target" wavefunction. Fig. 3 shows the evolution 
of this "overlap" , in terms of the scalar product squared. A relatively fast (tsweep = 7500 time units) and a slower 
sweep {tsweep = 15000) are shown. While the final overlap is about 0.9 for the slower sweep, it only reaches 0.7 for 
the faster sweep. Table II shows this final overlap for a choice of parameters. 

From Eq[2]one can form the adiabaticity parameter ^lunnei^ sweep / f-, which is shown in the next-to-last column of the 
Table. This parameter, which compares the splitting at "crossing" with the "velocity" e/tsweep of passing through the 
crossing, will be recognized as that which arises in the Landau-Zener theory As would be expected, one observes 
that values of the parameter of order one characterize the transition from adiabaticity to non-adiabaticity, and that 
similar values of the parameter give similar results. However, there appear to be some small deviations, (see the 0.90 
values), showing the usefulness of detailed numerical calculations. A value of about 4 or greater for the adiabaticity 
parameter appears necessary to achieve a 90% success probability, while a value near one gives the 50% point. With 
the typical time unit of 6.3 x 10~^^s a sweep of thousands of time units corresponds to some nanoseconds. Table II 
exhibits the great sensitivity of totunnei to /3, a reflection of the exponential nature of the tunnel splitting. 

In the simulations the sweep is carried out by simply letting 0"^^* vary at a constant rate as it goes from some initial 
value to the same value with the opposite sign; that is, the sweep is simply linear with a constant "velocity" . In more 
refined versions of the sweep procedure, it is conceivable that special waveforms could be used so as to pass through 
the dangerous vicinity of 0'^^'=G slowly while the overall sweep is very fast. However, it should be kept in mind that 
the decoherence is also most effective in the vicinity of "level crossing" (see below). All calculations presented here 
are performed with simple linear sweeps. 



DECOHERENCE 

We shall present a simple model for the decoherence suitable for numerical simulation and then examine its effects 
in various contexts. We first apply the model in the "spin 1/2 picture" where it can be handled analytically and then 
compare with numerical calculations with the full Hamiltonian. 




8 



13 




initial 0'="^* 


e 


^tunnel 


t sweep 


2 , —1 


final overlap 


1.14 


16.3 


0.002 


0.055 


0.027 


1 000 


13 


1.0 


1.14 






0.055 


0.027 


300 


4.0 


0.90 


1.18 






0.058 


0.0065 


10 000 


7.3 


0.98 


1.18 






0.058 


0.0065 


4 000 


2.9 


0.78 


1.19 






0.060 


0.0044 


15 000 


4.0 


0.90 


1.19 






0.060 


0.0044 


7 500 


2.4 


0.70 


1.19 






0.060 


0.0044 


4 200 


1.4 


0.50 


1.21 






0.063 


0.0019 


80 000 


4.6 


0.90 


1.21 






0.063 


0.0019 


40 000 


2.3 


0.68 


1.29 






0.075 


0.000045 


120 000 


0.0032 


0.016 



TABLE II: Tunnel splitting uotunnei as a function of /3 and its effect on adiabaticity. One observes a rapid change in ujtunnei 
with 13. The resulting effect on the "success" of an adiabatic inversion is measured by the "final overlap" as in Fig. 3. Similar 
values of the adiabaticity parameter give similar results for the overlap. With our "typical parameters" one thousand time 
units is 6.3 ns. 

Introduction 

To discuss decoherence it is necessary to introduce the density matrix p For our present purposes we think of 
/9 as a matrix arising from an average over different wavefunctions; 

p = ^^^ = -Y^^-^-\ (10) 

a=l 

where each ip'^ is a column vector and so p is a matrix. Our model will simulate the decoherence as a noise in the 
Hamiltonian, producing different wavefunctions in the evolution from an initial time to a final time. We then average 
over these wavefunctions according to Eq[TO]to obtain the density matrix. 

We recall the description of decoherence for the two-state system when thermal over-the-barrier transitions are 
neglected The decoherence is characterized by the parameter D arising in the effective "Bloch equation" giving 
the evolution of the density matrix. That is, the 2x2 density matrix is written in terms of the Pauli matrices a 

p=i(l + P-a), (11) 

where the information about the state of the system is in the "polarization vector" P, the average value of the abstract 
"spin", P = Tr[pcj]. 
The time dependence of P is given by the equation 

P = PxV-i:iPr. (12) 

where "T" stands for "transverse" , and means the components of P perpendicular to the direction in the abstract 
space chosen by the external perturbations. In the following we take the outside perturbations to be along the abstract 
z-axis associated with the basis states defined earlier. Pt then refers to the x,y components and represents the degree 
of quantum phase coherence between the two basis states u± . The shrinking of Pt induced by D in Eq [12] signifies 
a loss of phase coherence between the basis states. The parameter D thus characterizes the decoherence rate of the 
system. We shall also use tdec — 1/D to refer to the decoherence time. 

An important conclusion one can draw from Eq[T2]is that the major contribution to the decoherence occurs at "level 
crossing" . In a sweep passing through a "crossing" the P vector swings from "up" to "down" so that at "crossing" it 
is purely horizontal, with a large Pt- If we take the scalar product with P in Eq[T2l 

i^p2 = P • P = -DP • Pt = --DPt^ ■ (13) 

Since the departure of |P| from 1 measures the loss of coherence, the equation shows that the greatest "shrinkage", 
i.e. loss of coherence, occurs when P is transverse, at level crossing. 
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Random field in the full Hamiltonian 



To simulate the decoherence numerically we adopt a random field approach, where we suppose a small random 
time-dependent perturbation present in the Hamiltonian. This gives difi'erent realizations a of the Hamiltonian iJ". 
Starting with a given initial wavefunction, we evolve it with H°'(t) to yield a final "0°. Multiple repetitions of this 
procedure, with an average over the different realizations "0" according to EqllO( gives the density matrix originating 
from the initial pure state. This procedure can be implemented on the computer in a straightforward way. 

We shall model the perturbations due the external environment as a kind of fiux noise, with a random noise M 
added to (j)^^^ such that 

0'=^* f/.'^^* +A/'"(<) . (14) 

The external flux is of course not the only possible source of noise, in principle one may envision fluctuations of any 
of the parameters in the Hamiltonian. However the flux noise may well be a good representation of the decoherence 
and may also correspond to the main source of external noise in actual experiments. Naturally, many of our general 
conclusions will remain valid regardless of the specific origin of the noise/decoherence . 

Neglecting TV^ effects, Eq[llimplies - 0'=^*)2 i(0 - 0'=^*)2 - (j)M°-{t), so that the Hamiltonian becomes 

H' = ^^ + M\U - 0'^"*)'] + - Vocj^M^it) . (15) 

To interpret the additional term we use the identification Eq|4l relating to the of the spin picture: 



- Vo<j>M^it) « Fo0,AA''(t)a, = a^B^it) , (16) 

where we call the quantity — Vb^cA/" the random field B. Thus in the spin picture the noise term has the interpretation 
of an additional random "magnetic field" B applied to the z- component of the spin. We take this random field to 
have average value B — 0. 



Random field in the spin picture 

We begin with an analysis within the spin 1/2 picture, which can be handled analytically and which provides an 
orientation for the full simulation. The Hamiltonian of the spin picture is now 



H'' = ^(7 -V + B^a, . (17) 

We first try to establish the connection between the decoherence parameter D in Eg [T^ and the statistical properties 
of B. Writing a wavefunction in terms of "up" and "down" basis states u± as tjj — au^ + the contribution to 
the density matrix from one instance in Eg 1101 is 

p-(:f. . (18) 



a* (3 (3(3 

where aa* + (3(3* = 1 from normalization. Taking the average and comparing with Eq llH one has 



= Re a(3* Py = Im a/3* (19) 

We now assume that the frequencies associated with the noise are much above those connected with the slow coherent 
rotations induced by V, so that V maybe thought of as effectively "turned off" for the calculation of D. In this case 
the solution of the Schroedinger equation is simple. If at i = we have the state ip = aou+ + (3oU-, then at time t 

^"(t) = ao elo ^^'''u+ + (3o e" lo ^^"^'u^ , (20) 
Putting q;o and (3^ real so Pt is initially along the x-axis one has according to Eg [TOl 
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= (ao/3o)Re e'' L (21) 
Assuming the random field is of small amplitude permits the expansion 



= («o/3o)Re e'' Jo ^"^ « (ao/3o)(l - 2 ( / Bdi)') = («o/?o)(l - 2 / / dt" dt' B{t")B{t') ) , (22) 

Jo Jo Jo 



Jo 



where the linear term vanishes with B = 0. We thus have to do with the autocorrelation function B{t")B{t') — 
B'^{t")B'^{t') where i?° is a particular realization and N the number of realizations. 

As is known in the theory of random noise [l^ a correlation function of this type under stationary conditions 
behaves such that /g /g dt" dt' B{t")B(t') = t B{t)B{0)dt = 2tJ^ B{t)B{0)dt. 

Now with V = and using, as follows from Eq[T2 that — Pxoe^^* ~ (ao/3o)(l ~ Dt), comparison with Eg 
yields the identification 

D^-2[ [ dt"dt'B{t")B{t') = 4 / B{t)B{Q)dt . (23) 
i Jo Jo Jo 

Thus in the "spin 1/2 picture" D is given by the integral of the autocorrelation function of the random field. Or, 
regarding azB"" as a perturbing energy 6E, one can also say D = 4 dE{t)6E{0)dt. 

As is evident from this short derivation, or from the assumptions used in the original derivation [8] of Eq 1121 it is 
assumed that the frequencies of the random perturbations entering into D are high compared to the slow coherent 
motions induced by V. In the thermal context, where one anticipates the perturbing frequencies to be on the order 
of the temperature, this means that the energy splittings induced by V are assumed small compared (k=l units) to 
the temperature. For this reason it is consistent that in Eq[12]the density matrix relaxes to the identity and not to 
the form that would be given by a Boltzmann factor. Similarly, low frequency instrumental noise in the laboratory 
would be better treated as an additional contribution to V rather than being incorporated into D. In our simulations 
we always treat the noise as being of high frequency in this sense. 

Modeling of the noise 

A simple noise model suitable for numerical simulation is 



Mit) - 7^it)A , (24) 

where rj is a random sign rj = ±1, and A a positive magnitude. Let St be a certain small time interval during which rj 
is constant and let the probability that there is a sign switch for the next interval be psw (and to remain unchanged 
1 — Psw)- This procedure, in the limit of small Psw and St, leads to an exponential distribution for the noise pulse 
lengths and an autocorrelation function 

7]{t)rj{0) = e-2P=™*/'5* = e-2«* (25) 

where a = Psw/St. Introducing the noise power to characterize the frequency content of the noise signal [lo| . such an 
autocorrelation function leads to the noise power spectrum dt cosLut 77(^)77(0) = ZJ^^^^ with Wc = 2a = 2psw/St 
a cutoff frequency. This spectrum is roughly constant ( "white noise" ) up to the cutoff at about tOc and then falls off 
as at higher frequencies. We attain the highest cutoff, the closest approximation to infinite frequency white noise, 
by choosing the switching time St as small as possible i.e., equal to the program step and with the switch probability 
Psw = 1/2. In this case ujc = 1/ program step. 

For Eq[23]we need the time integral of the autocorrelation function, which is 

POO 

/ M{t)M{0)dt = AVwc , (26) 
Jo 

Using Eq [23l and recalling the definition of B in Eq[T6l one obtains, finally 



i^ = 4(Vb0e)' — 



(27) 



11 



The 1/lOc behavior originates in the fact that one is essentially performing a random walk in the phase of the 
wavefunction, and the step length of this random walk is ~ A (5i ~ A/wc- In the next section we check this prediction 
from the spin 1 /2 picture against simulations with the full Hamiltonian. 

Tests of the noise simulation with the full Hamiltonian 

We now proceed numerically, using the SQUID Hamiltonian Eq[T5] . An initial state is evolved with Eq [15] with a 
given noise realization A/"" . Repeating this procedure many times a final density matrix is obtained by averaging over 
different realizations according to EqlTOl If the two-state approximation is good, we should observe a decoherence 
rate in agreement with Eq l27l 

The numerical calculations were again carried out by repeated iterations of (1 + iHAt). As few as 30 samples of 
the random field often give good (10%) results. A time step Af « 1 was usually found sufficient. For higher accuracy 
more samples and smaller time steps can be used. 

It should be noted that the density matrix arising here from Eq [10] has a slightly different significance than in the 
"spin 1/2 picture". Since we now work with full wavefunctions 4' (4') in the "position coordinate" (f), all eigenstates of 
the SQUID Hamiltonian are potentially present in the wavefunction and in the density matrix p((/)',<^). Therefore, to 
compare with the results of the previous section we must specify some basis of wavefunctions and evaluate the matrix 
elements of p{(l)',(f>) in that basis. We shall use cither the two lowest energy eigenstates, or the "up", "down" basis 
corresponding to the eigenstates of Cz. 
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FIG. 4: Effect of decoherence, as simulated by random noise in the full Hamiltonian Eq 1151 for the symmetric case with 
^ext_Q SQUID parameters were /3 = 1.19,^ = 16.3. The initial state is the lowest energy eigenstate. The pure state at t=0 
is converted to the maximally mixed state. The quantity plotted is pn in the basis of energy eigenstates, which for (j)'^^'^—0 
corresponds to ^(1 + Px)- The data is fit to |(1 + e~^*) yielding D = 0.0017, in good agreement with the prediction D=0.0018 
from Eq[27l Noise parameters were A = 0.00032, = 0.05. With the "typical" time unit {L/pH C/pF)^^^ x 1.0 x 10~^^ 
seconds^ 6.3 x 10~'^^seconds, the decoherence time 1/D = 560 corresponds to 3.5 ns. 



Damping of Pt 

In general the evolution of the density matrix will reflect the simultaneous effects of the internal Hamiltonian and 
the external interactions. In Eq 1121 D gives a damping of the transverse components of P. Thus the simplest way to 
see the effects of the decoherence alone is to start from an eigenstate of the symmetric, (/)'^^*=0, Hamiltonian where 
V is parallel to P and purely transverse. In the "spin 1/2 picture" P would start in the x direction with value 1 and 
decay exponentially with decoherence time 1/D to the value 0.5. 

For the numerical simulation of this situation with the full SQUID Hamiltonian we begin with the lowest energy 
eigenstate, obtain the density matrix at a certain time, and take its expectation value with respect to the starting 
wavefunction, which quantity we call pu. Fig. 4 shows the results using the noise parameters A — 0.00032, ojc — 0.05 
for Af. An exponential fall-off is observed, and a fit gives D = 0.0017. From Eq [^7] with (f>c ~ 0.90 we predict 
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D = 0.0018. This is in good agreement with the fit. We note that this value of D is not large in comparison to 
the energy sphtting 14 = 0.0044, so that in Eq [12] both terms on the right are of the same order of magnitude. A 
prediction of the curve from Eq[T3 using this D and is essentially indistinguishable from the fit curve. Excitations 
to higher states were allowed by the program ("second cut" =4), but none are evident in that the decay of the curve 
is to 0.5 and not a smaller value. Runs with two or four states for the "second cut" showed no significant differences. 



Damped Oscillations 

To exhibit the characteristic two-state oscillations, Fig. 5 gives the results of a simulation run with the same 
parameters, but now with the initial state chosen to be an eigenstate of <Tz, i.e localized in one of the potential wells. 
The eigenstate of is found by constructing the 2x2 matrix for (p in the two first energy eigenstates and then finding 
the eigenvectors of this matrix. 

We anticipate that the wavefunction will oscillate back and forth with a damping governed by D, as is indeed seen in 
Fig. 5. The quantity plotted is the density matrix element with respect to the starting state, which has the two-state 
interpretation pn = (1/2)(1 -I- Pz{t)). Again we find no perceptible difference between runs using a two-state and a 
four state basis for the time evolution, indicating little excitation of higher states with these parameters. 

1 H , , , , , , , , , 1 

\ 

0.9 -'i 



0.8 - 




1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 
Time 



FIG. 5; Same conditions as Fig. 4 but with the initial state chosen as an eigenstate of tj^. Since this is not an energy eigenstate 
oscillations occur, which are then damped by the decoherence. The vertical axis corresponds to ^(1 -I- Pz). The curve is to 
guide the eye. 



Turing-Zeno- Watched Pot Effect 

We briefly look at strong damping, which is the limit 

— ^»1. (28) 

^tunnel 

Note that in dividing Eq [T^] by V ^ one obtains an equation containing only the scaled time variable tV and the 
parameter D/V, so that when using a time variable appropriately scaled to the oscillation frequency, D/V is the 
only parameter in Eq[T21 With large D/V we enter the regime of the "Turing-Zeno- Watched Pot Effect" where one 
expects that the damping D inhibits the tunneling strongly; if the system is in one of the potential wells it tends 
to remain there. We approach the "classical" situation where quantum mechanical linear combinations cease to exist. 
Fig. 6 presents some simulations of this situation. The conditions are as in Fig. 5, but with increasing A. Runs with 
other oscillation frequencies, that is with different /3, yield the same results when the time is appropriately rescaled. 
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FIG. 6: "Watched Pot" effect. Same conditions as in Fig. 5 but with various degrees of strong decoherence D/uJtunnei » 1. 
The "freezing" of the evolution is observed as predicted Q| from Eq[T2]with p(l, 1) « 1 — \{<^tunnei/ D)t. The time interval 
shown would correspond to about seven undamped oscillations. The different values of D were produced by varying A from 
0.01 to 0.002 and calculated from Ea[27l 

Noise Frequency 

According to Eq[571 changing and lUc such that their ratio remains constant should leave D the same. However, 
the higher frequency components of the noise spectrum may have an independent effect through the excitation of 
higher states. Excitations beyond the lowest states would be undesirable, as it represents a non-unitary evolution 
among the lowest states, with some probability going to higher states. 

On the other hand, if the frequency of the noise is much less than the frequency corresponding to the distance to 
the next set of levels, the adiabatic theorem tells us that the states remain in the lower set and that such non-unitary 
effects are suppressed. Since our principal level splittings are generally substantial fractions of unity, there should be 
little excitation of higher states for noise frequencies << 1- In a plot of the type Fig. 4, excitation of higher states 
is manifested by the decay of the density matrix element to a value less than 1/2, showing population beyond the 
first two states. 

To examine these points, we show in Table III a series of runs at constant A^/wc, but different oJc- The resulting 
D as determined from a fit as in Fig. 4 and the final p(l, 1) are shown. The approximate constancy of D seems to be 
well verified. For << 1 the decay is indeed to 1/2, but for larger values there are noticeable departures. In another 
set of runs with A reduced by a factor of two this effect is almost entirely absent. With our typical time unit these 
values of D corresponds to a decoherence time of (1/-D) x 6.3 x 10"^^s = 3.5ns and 14 ns. It should be noted that 
our noise spectrum has a relatively strong high frequency tail ~ 1 /uj"^ as opposed to the exponential cutoff one would 
expect for a purely thermal background. 

The actual value of D is of course our great unknown, and we treat it as simply a phenomenological parameter. 
For orientation we can keep in mind the estimate D = T / [e^R) which we have used previously; or measurements by 
direct observation of the damped oscillations ^HJ on a similar SQUID system. These find a decoherence time « 20 
ns at 25mK. This corresponds to the reasonable value R = 260 kfl for the effective resistance. Our sample values of 
tdec = 1/ D ^ 10^ in the dimensionless units, with our typical time unit of 6.3 x 10~^^sec, would correspond to some 
nanoseconds or tens of nanoseconds for tdec- 

DECOHERENCE IN THE NOT GATE 

Having checked that our noise/decoherence simulation has reasonable features, we now turn to some applications. 
The one variable or one-qubit problem is the simplest situation. When (j)'^^*' is swept (from a relatively large value in 
the sense Vz » Vx) to its opposite value, interchanging "up" and "down", it represents the logic gate NOT. The 
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A 


D, fit 


D, EqM 


final p(l, 1) 


2.0 


0.0020 


0.0017 


0.0018 


0.35 


0.50 


0.0010 


0.0017 




0.40 


0.11 


0.00045 


0.0018 




0.48 


0.050 


0.00032 


0.0018 




0.50 


0.025 


0.00023 


0.0019 




0.50 



TABLE III: Values of D and final p(l, 1) for runs as in Fig. 4, with varying loc but constant A'^/loc- The prediction of an 
approximately constant D, as well as the value for D, are in agreement with Eg 1271 The relaxation to a value below 1/2 for 
the larger A and luc values indicates excitation of the third and fourth levels. 
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initial 0'=^' 




A 


t sweep 


final overlap 


1.19 


16.3 


0.0020 


0.042 


0.000042 


30 000 


0.95 












60 000 


0.96 












80 000 


0.90 












150 000 


0.85 












800 000 


0.60 












1 000 000 


0.58 












2 000 000 


0.52 



TABLE IV: Effect of decoherence on the "success" of adiabatic inversion sweeps. Noise parameters are chosen to give a 
decoherence time tdec — 1/D = 29000 time units. "Final overlap " is the average of the overlap squared -^Sa] < V'^IV'/ > 1'^ 
or equivalently < V'/IpI^/ >• ^.s the sweep is made slower, so that decoherence has time to take effect at level crossing, 
p approaches the totally incoherent value of 1/2. Errors on "Final overlap " are on the order of a few percent. For the 
no-decoherence situation see Fig. 3 and Table II. 

understanding of the effects of decoherence here is important both for finding the regime of operation of the logic gate 
and for our proposal fll] for measuring D via the success or failure of an adiabatic inversion. 

In some applications of the adiabatic idea in quantum computing a special role is assigned to the ground state, as 
in the search for the minimum of a complicated functional [l3|. However our simple gates are not of this type and 
both states of the qubit, either the ground state or the first excited state, are on an equal footing. Thus in the NOT 
operation, for example, it is equally important that 1 — > or ^ 1, and which of these is represented by the ground 
state is of no particular significance. 

As in the discussion concerning adiabaticity, we measure the "success" of a sweep by the value of "final overlap" as 
in Fig. 3, where now there is an average over realizations of the random noise. Table IV shows the effects of increasing 
sweep time with fixed noise parameters. The sweeps are sufficiently long, according to Table II, that non-adiabaticity 
should be unimportant. 

Interestingly, although the noise parameters have been chosen to give a decoherence time of 29 000 time units. 
Table IV shows decoherence not fully setting in until significantly longer sweep times. This may be understood in 
terms of our remarks in connection with Eq[T3]that decoherence has most of its effect during level crossing. If one 
attempts to estimate the time spent during the sweep in the vicinity of the "crossing" with the given conditions, it 
is on the order of 10%. This suggests that the relevant time for the onset of significant decoherence is not 29 000 
time units but rather more like 290 000, in agreement with the behavior in Table IV. This argument is supported by 
examination of Fig. 3 or Fig. 13, where one sees that the actual switching of a state takes place in a small fraction of 
the total sweep time. Hence for our sorts of logic gates, the speed needed to avoid decoherence may be less stringent 
than one might simply infer from comparing the sweep time to the decoherence time. 

To illustrate the effects of non-adiabaticity and decoherence together, we show in Fig. 7 a series of simulations for 
varying sweep times and A's. There appears to be a broad range around twenty thousand time units where neither 
effect is drastic. The times refer to our dimensionless units, so for our typical time unit 6.3 x 10~"'^^s, the two values of 
D in the plot would correspond to 1/D — tdec =50 ns and 180 ns. If one were to consider SQUID parameters giving 
a larger ujtunnei and so a shorter adiabatic time, the favorable region could be widened considerably to the left, i.e. 
to shorter times. (See the discussion below "Smaller /3") 
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FIG. 7: Adiabatic and decoherence effects on the "success" of a sweep (—0.002 0.002) in (/!>'^^' , for the one dimensional 
system (NOT) as the sweep time is increased. The adiabatic time for 50% "success" is indicated by an arrow. As the 
sweep time is lengthened one observes the onset of adiabaticity, followed by decreasing "success" as the decoherence takes 
effect. Squares and dots: noise parameters adjusted to give D = 1.3 x 10~^. Squares: lOc — 20mK and A — 5.0 x 10~^. 
Dots: Uc = 5Qmii'; A = 7.9 x 10"^ Triangles: noise parameters adjusted to give D = 3.5 x 10 using Uc = 50mK and 
A = 4.2 X 10"^ SQUID parameters are (3 = 1.19, C = 16.3. 



TWO VARIABLE SYSTEM 



We now turn to the two variable or two-qubit system. When two one variable systems are weakly coupled we arrive 
at a two variable system with four low lying states. For the SQUID these would be the four possible states arising 
from the current circulating clockwise or counterclockwise in two SQUIDs, and as explained in ref Q this can be so 
arranged that the result of the sweep of one of the (Z)*^^* depends on the state of the other SQUID, thus providing the 
conditions for a CNOT gate We recall 0] the Hamiltonian for this problem 



-1 -1 

^^^^ + ^^ + ^ ^ (29) 
2/ii d(l)i 2^2 0(/)2 



with V 

V = - 0r*)' + h{cl^2 - (^Tf - 2^12(02 - - (I^T)] + /3lC0S(/.i + ^2003(^2} • (30) 

In place of /i = Vq for the one variable case, one now has ~ Vq. Analogously to the single variable case, 

the factor Eq which converts the energy of the dimensionless Hamiltonian to physical energy involves the inductance 
and capacitance of the two SQUIDs 2], namely Eq = 1/^LC, where — i^l^^^^-i - and C = \fC\C2- The small I 
are the inductances refered to L. In Fig. 8 we reproduce a contour plot for the potential showing its four minima, 
corresponding to the four states of the 2-qubit system. 

As in the one dimensional case, the numerical calculations use a large harmonic oscillator basis, now in two variables. 
For most runs a basis of one or two thousand states was used. Even when using the "second cut" reduced basis to 
four eigenstates, time dependent runs with many samples could take several minutes on fast PC's. It seems that 
extensions to more than two or three variables would need new computational methods. 
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FIG. 8: Potential contours for Eg 1301 showing the four minima of the 2-qubit system. 

A basis for the four states is provided by "spin up(l)spin down(2)=?i_|_(l)u_(2) and so forth: 

u_(l)ii_(2) u+(l)ii_(2) w_(l)u+(2) u+(l)u+(2) . 



(31) 



These four states correspond to definite senses for the currents in the SQUIDs. 

On the other hand if the parameter I12 giving the flux coupling between the SQUID is very small, the situation 
reduces to two independent single variable systems. Thus for the analysis of weak coupling with small /12 it is 
convenient to introduce the eigenstates for any value of the individual (j)'^^* for each SQUID: 



^-(l)^^-(2) 



v+{l)v-{2) 



«-(lK(2) 



v+{l)v+{2) 



(32) 



This is the independent SQUID basis, where each SQUID can have its own Hamiltonian, according to its t/)'^^*. 

The V are those discussed in connection with Table I: u+ — (cos| u+ + sm| m_) and v- = (cos| u_ — sot| u+) . If 
the applied (ff^^aie different, then we have different angles 9i and 02 in these relations. At I12 = the u(l)w(2)'s are 
the eigenstates of the complete system, and with I12 7^ there will be a mixing among them. 

COUPLED HAMILTONIAN IN THE SPIN PICTURE 

In the spin picture for Eq [30] there are now two "spin 1/2" objects, interacting through I12 and subject to the 
external fields 0^^* , 02^* . By the arguments used for Eq Hj we make the identifications 



01 ^ 0c(l)fT,(l) 

and the effective spin Hamiltonian is 

i a{l) • Vi i a{2) ■ V2 - ^0^12 ( 0c(l)a.(l) - ^f* 



Cj)2 0c(2)(T.(2) , 



0e(2)a,(2)-0^ 



(33) 



(34) 



For small I12 the components of the V are determined as in Eq[Sl namely Vz ~ 2Vo I (j)c4>'^^* while 14 is found via 
Eq[71 As before, (pc is approximately the of a state localized in one of the potential wells. The effective interaction 
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FIG. 9: Left: Square of the wavefunction for state 1 localized in the lower left potential well of Fig. 8. Right: at "level crossing" 
as the wavefunction moves to the lower right potential well. Potential parameters were h — h ~ 1, li2 = 0.005, /3i = /32 = 1.19 
and Vb ~ 16.3. 

Hamiltonian between the two devices is then 

~ -Voh2M2 ~ -l/o?i20c(l)0c(2)a,(l)a,(2) . (35) 

This operator induces level shifts of the four u base states without mixing them. However when the tunneling is 
introduced as a perturbation, non-trivial combinations of the four states can arise. 

Low-lying level patterns 

We begin with weakly interacting SQUIDs, I12 << h.h- The energy level pattern anticipated for the lowest levels 
may be understood by beginning with the totally decoupled system of just two independent devices, I12 — 0. If we take 
both devices with about the same parameters, for example, we have first the ground state where both SQUIDs are in 
their lowest state, the first state of Eq[32l Then there are two approximately degenerate states with one SQUID in the 
first excited state and the other in its ground state, the second and third states of Eq[32l Finally the fourth state has 
both SQUIDs in the first excited states, the last wavefunction of Eq[32l If each device is in the configuration where 
the splitting of its first two states is small (small 14), as in our discussions above, then the splitting from the ground 
state to the degenerate pair is equal to the splitting from the pair to the fourth state. Furthermore the splitting to 
the fifth state should be substantially greater since it involves an excitation of the principal quantum number. 

Turning on a small we show an example from numerical calculations in Table V. With I12 = 10~^ these reveal 
the expected pattern. One sees that the fifth state is well separated from the lower ones, again supporting the use of 
the picture of an approximately isolated Hilbert space, as for the first two states in the single variable case. 

In the example of identical parameters for both SQUIDs, with second and third levels degenerate for I12 — 0, 
we may use degenerate perturbation theory to find the splitting induced by turning on I12. The matrix element is 
< v+{\)v-{2)\Hint\v-{l)vj^{2) >= Voli2{sin9)'^ , where 9 is the angle for going from the u to the v as discussed in 
connection with Table I. The splitting of levels 2 and 3 is then 2Voli2{sin9)'^ (j)'^. With our typical parameters and 
Z12 = 1 X 10~^ and (j)'^^^—0 and so sin9 = 1, the formula, using (pc ~ 0.88, gives 25 x 10^^. The numerical calculation, 
as shown in Table V, gives 27 x 10~^. 

Observation of this small splitting would be quite amusing with regard to the question of coherence between 



macroscopic objects. While work with the "Cooper-pair box" [1J| has seen effects involving interference between the 
states of two qubits, the states involved differ only microscopically, namely by a Cooper-pair. Here, with the SQUID, 
the states concerned differ by the circulation direction of a macroscopic number of electrons, and so seem to pose 
the macroscopic coherence question more dramatically. The energy eigenstates resulting from the diagonalization 
to obtain the small splitting are (l/-y2)(w_(l)w+(2) ± z;_|_(l)i)_(2)). The splitting, one might say, results from the 
relative quantum phase ± involving the two SQUIDs, that is, between two macroscopic objects. This is yet a step 
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FIG. 10; Level behavior for the first four states of two identical SQUIDs with coupling Z12. In the first panel one notes the 
pattern of a ground state, an almost degenerate pair and single state, as would be expected for two independent SQUIDs. This 
pattern changes as the interaction energy between the SQUIDs ~ V0I12 becomes on the order of the splitting cotunnei- 



further than the effects for one SQUID, where one is sensitive to the phase between different macroscopic current 
directions, but in one device. Since this splitting is very small, however, line broadening due to the noise effects may 
be comparable to the splitting. Values oi D ^ 10^* — 10~^ for the tens of mK region could be in the same range as 
the splitting for I12 ^ 10^^. 

The degenerate perturbation theory formula only applies when we have the picture of two approximately degenerate 
levels well separated from the others, as in Table V. The picture changes rapidly as the mutual coupling I12 is increased. 
In Fig. 10 we show the energy level patterns for the first four levels for different I12 as functions of a common applied 
fj^ext^ Again, the two SQUIDs are taken to be identical, with our standard parameters. One observes that after about 
li2 ~ 10~* the uncoupled SQUID pattern no longer holds, since now the interaction energy ~ V0I12 is on the same 
order as the original splittings, oJtunnei = 0.0044. A characteristic change of behavior occurs when I12 «(/)'^^'. This 
may be understood as the value of I12 where the flux contributed from the other SQUID becomes comparable to the 
applied (f)'^^^ . (see Eq 5 of rcf 3] ) . 

Another limit which is not difficult to analyze is that of small 9 or relatively large (jf^*{see Table I). In this limit 
the u's are the eigenstates and the interaction Eq [35] simply gives additive contributions to the energies. We may 
even consider different parameters and (ff^^ for the SQUIDs. According to Eq [34] this results in a ground state 
with both spins "down" and energy i(— 1/2(1) — 14(2) — 5) and an upper state with both spins "up" and energy 
\{+V^{l) + K(2) - 5). There are two middle states with energy \{+V^{l) - K(2) + 5) and \{-V^{l) + V^{2) + 5); 
5 is the contribution from Eq[3ni <^ = ^'i2'/'c(l)'/'c(2). In the case of identical SQUIDs these middle states form a 
degenerate pair, as one sees for the larger (/)^^*. With small deviations of the V from the z-direction there will again 
be a splitting which can be found from perturbation theory. 



CNOT CONFIGURATIONS 



Our design for a CNOT operation consists of an adiabatic sweep from a point with relatively large (/)'^^* in the 
(05^^*, 01^*) plane to another such point so that there is a definite mapping among the four u(l)u(2) states. The 
four possible configurations of currents clockwise and counterclockwise in the two SQUIDs undergo a certain definite, 
reversible, rearrangement. This rearrangement is chosen in accordance with the definition of CNOT: one pair of 
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TABLE V: Splittings at (jf^^—Q among the first five levels for the first panel of Fig. 10. The energy difference AS is with 
respect to the previous level, thus the first excited state is 0.0043 units above the ground state. The fifth state is distinctly 
separated from the first four. The first and third splitting should be approximately that of a single SQUID, which with (/)'^^'=0 
is 0.0044. The small splitting between the second and third state is in agreement with the prediction of degenarate perturbation 
theory, 2^0^12 Re- 



states remains unchanged (control bit is zero), while the other pair reverses (control bit is one). This is accomplished 
by having the energy eigenstates (1,2,3,4), each one concentrated in a different potential well, move from one set of 
locations to another Since for relatively large (^"^^'each well represents a distinct state of the SQUID currents, one 
physical configuration of the two qubits is mapped to another. These rearrangements are represented by "tableaux" 

like tj] showing the localization of the first four energy eigenstates in the potential wells of Fig. 8. 



^1 2^ 

In a first step it is only necessary to identify those locations of the (0^^*, 02^*) plane with different tableaux such 
that a sweep from one to the other leads to the desired rearrangement. An example is^^ 2)^(^1 2^' "^^^ 

lower row may be identified with control bit zero, and the upper row with control bit one. The adiabatic theorem 
then guarantees that a "sufficiently slow" sweep will preserve the occupation of the levels, 1 — > 1, 2 ^ 2 and so forth. 
Once having so identified the desired sweep, the adiabaticity may be examined afterwards. 

Not all points of the {(j)'^* , (f)'^*) plane are suitable starting or ending points for a sweep since we require "well 
defined" wavefunctions. The wavefunctions should be A) well localized in one potential well so that the SQUID is 
in a definite flux state, and B) all first four states should be localized in different wells. Fig. 9 shows the appearance 
of a "good wavefunction " with state 1 in the lower left potential well of Fig. 8. In the searches a wavefunction was 
considered as "well localized" when the distance between centers of the wavefunctions was more than 2.5 times the 
spread of the wavefunction, as measured by the dispersion. 

In Fig. 11 we show the results of such a search. The associated tableaux are indicated for each region. The black 
areas are those of "bad wavefunctions " . The general range of "good" corresponds to the observations for one SQUID 
in Table I, where the spin picture was valid up to about (p'^^^Ki 0.01. 

In principle a CNOT may be accomplished by a sweep between any two adjacent regions where one row (or column) 
interchanges and the other does not. As was explained in ref. [s*], the switching values between tableaux (dark lines) 
may be understood as a "level crossing" occurring when the flux from the control SQUID onto the target SQUID is 
equal and opposite to the applied flux onto the target SQUID. That is, when (j)'^^*= l^/l- At this point the total flux 
on the target SQUID is approximately zero and so is at a level crossing. This argument allows one to understand the 
pattern of dark lines. It will be noted that in crossing a dark line only adjacent energy levels switch, as expected for 
a "crossing" . The vertices are singular points involving more than two levels. 

For a smaller I12 the map will then be similar but with smaller regions between the vertical and horizontal dark 
lines. In the central box of the plot, although there are "good" wavefunctions, the applied flux is apparently too 
low for the "immobilization" argument of ref. [3:| to work and there are no CNOT rearrangements within the box. 
However CNOT sweeps exist connected to peripheral regions. These take place in the vertical or horizontal direction 
with the flux on one qubit (control bit) relatively large and constant while the other flux (target bit) is varied. 

For such vertical or horizontal sweeps, involving changing one flux, only one SQUID inverts, according to which one 
undergoes the sweep. In diagonal sweeps involving both fluxes on the other hand, states transfer across the diagonal 
in the tableaux, indicating changes in both SQUIDs. This explains why the diagonal black lines are very narrow, 
since such "double flips" involve two tunnclings with a corresponding small mixing energy. 

Although the map of Fig. 11 was obtained by an examination of the individual wavefunctions, a reasonable idea 
may be had by simply inspecting the potential landscape. Since the ordering of the energy levels will usually follow 
the ordering of the minima of the potential wells, we find that the ordering of the minima usually give the correct 
tableaux. Thus a suitable sweep is frequently simply one where the potential minima rearrange in the desired manner. 
Alternatively, all the transition lines on the map may be found from the set of linear relations arising from the "level 
crossing" conditions (coefficients of the az) while treating the tunneling (a^) terms as a perturbation. 
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FIG. 11: Map of the regions of {cpT i^i )plane with well defined wavefunctions suitable for performing CNOT. The numbers 
in a tableaux show the location of the first four energy eigenstates on the potential landscape of Fig. 8. These show that a 
CNOT may be obtained by horizontal or vertical sweeps in the regions external to the central square. Potential parameters 
used were h = h = 1, h2 = 0.005, /3i = /32 = 1.19 and Vo = 16.3. 



Hilbert space completeness 



For the "spin 1/2 picture" to be a meaningful in a complex system with many levels it is necessary that the 
states selected constitute an approximately independent vector space. If, for example, the selected states were to 
mix with other states of the system as we carry out our gate operations, then the two states would not be a faithful 
representation of a two-state qubit, since more than two states would be involved. And similarly when we have two 
qubits or four states for CNOT, these four should act effectively as a separate space. 

The Hilbert space completeness test mentioned at the end of the section "Identification with a spin 1/2 system" 
proves to be quite interesting in this regard and may be a useful tool in studying higher dimensional systems. In this 
method, we take a point in the parameter space of the Hamiltonian and compare it with some reference point. If 
the Hilbert spaces for these different Hamiltonians are closely the same, there should be a high overlap of the lowest 
eigenstates for the different Hamiltonians. As we shall explain, a test can formulated which applies to any linear 
combination. For the one variable system the parameter space in question is the 0*^^* line, and with two variables the 
((^f^*, 02^*) plane. The precise point we choose as the reference is unimportant in regions where the overlaps are high; 
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FIG. 12: Region of the {4>T"^,4>T*) plane with high "completeness " using Eg 1361 One notes the close similarity to the region 
defined by Fig. 11, which was found by examination of individual wavefunctions. 

for the *) plane we take (^f *, 0^^*)=(O,O). 

Let the Vi be the lowest eigenstates at the reference point and v' a linear combination of the lowest eigenstates for 
some other parameter values. We search for the "worst case", the minimum value of E| < Vi\v' > p. The procedure 
involves constructing the matrix My =< Vi\vj >. For the one SQUID or two state system i and j run from one to 
two; in the two variable system from one to four. 

One now observes that the sum of the overlaps squared Sj < Vi\v' > p is the expectation of the matrix M'^M in 
the state \v' >. Thus the worst or lowest overlap for an arbitrary normalized linear combination in the v' space is 
given by the lowest eigenvalue of this matrix: 

least overlap = smallest eigenvalue M'^ M (36) 

Varying the chosen point in the plane, one maps out the region of a common Hilbert space where the 

worst overlap is close to one. Fig. 12 shows such a map, produced with high resolution. A grey scale was used for the 
points, from white for least overlap = 1 to black for least overlap = 0. 

It is quite interesting that the region of the (t/)^^', </)2^*) plane found by this method is essentially the same as that 
found by the detailed examination of individual wavefunctions as in Fig. 11. This may perhaps be understood by 
noting that if we have "bad wavefunctions " with two states in one potential well, this will have a poor overlap with 
the "good" situation where each wavefunction is in a different well. On the other hand "bad wavefunctions "in the 
sense that they are delocalized, as at a "crossing" , will still give a high "completeness" . The sharp transition, with 
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FIG. 13: Behavior of the states 1 and 4 in a CNOT sweep connecting two tableaux, where state 1 should remain in its original 
location while state 4 should reverse, as in the lower left of Fig. 11. No noise is applied. The vertical axis shows the overlap of 
the wavefunction with its intended final state. State 1 (upper curve) should not change its location and is seen to be stable. 
State 4 (lower curve) should reverse its location and is seen to go from zero to approximately complete overlap with the desired 
final state. State 2 and state 3 are not shown but behave similarly to states 1 and 4 respectively. Squid parameters are as in 
Fig. 11, and the sweep time was 20 000 units. It is to be observed that the switch itself takes place in a small fraction of the 
sweep time. 

little grey area, seems remarkable. Although the corresponding transition in the one dimensional case was rather 
abrupt, it seems more so here, presumably the effect is multiplicative. Both Figs. 11 and 12 contain 90 000 pixels, 
necessitating runs of many hours. 

Fig. 12 does not contain the narrow black bands of Fig. 11, where there are "bad wavefunctions " at "level crossings" . 
This is because these "bad wavefunctions " , although delocalized in and so not corresponding to a definite state of 
the SQUIDs, are still in the same Hilbert space as the localized states. An example of such a state is shown in the 
Right panel of Fig. 9. It should also be noted that the narrowness of the black bands in Fig. 11 is in accord with our 
earlier remarks concerning the small fraction of the time spent near "crossing" during a sweep. 



The extension of the previous noise treatment to the two variable problem is straightforward if we continue to 
represent the noise as fluctuations of the external fluxes 0*^^*. To each 0*^^* there will now be a noise term as in EgfTSl 
As in the discussion of Eq[34]we have Vz ~ 2Vb/(/)'^^*0c and so an additional term in Eg [34l 



Noise for the two variable system 



Fo{/i0c(l)fT.(l)AAi +/20c(2)a,(2)AA2} , 



(37) 



where again we keep only the linear term in the noise. 

There will then be a I? as in Eg [771 associated with each variable 



Di = 4 



(FoZi0c(l)Ai)' 

C^c(l) 



D2 



^(W20c(2)A2)^ 



(38) 



We can now carry out the simulations with the two noise signals imposed, one for each SQUID. 



CNOT SWEEPS 



As for the one-bit case, the "success" of a sweep can be measured by the overlap of an evolved wavefunction with 
the stationary eigenstate for the intended final state. For CNOT two wavefunctions, representing control bit =0, 
should not change their states, while those representing control bit =1 should move to new positions. We exemplify 
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FIG. 14; A CNOT sweep in the lower right ol the {(j>T\<j)T^) plane of Fig. 11, with (0, -0.01) (0.006, -0.01). The final 
overlap is plotted for state 1, which should change position. The coupling between the SQUIDs is I12 = 0.005. The sweep without 
noise (circles) shows the effect of non-adiabaticity for fast sweeps, while the addition of noise exhibit the noise/decoherence 
effects at longer times. The squares have noise parameters A — 0.000042, Uc — 0.042 and the stars A — 0.000079, Uc — 0.042. 
The noise is applied to both SQUIDs equally but independently. 



a successful CNOT sweep in Fig. 13 where we show the behavior of states 1 and 4 for a sweep between a tableaux 

In Fig. 14 we illustrate the effects of non-adiabaticity and decoherence for the 2)^(^2 l) '^^^^P '^^ 

right lower portion of Fig. 11, with (0,-0.01) (0.006,-0.01). The curves show the behavior of the overlap for 
state 1, which ideally would behave as the switching curve of Fig. 13, reaching close to complete overlap. The upper 
curve, without noise, shows the progression to adiabaticity as the sweeps become slower. The lower curves have noise 
applied, equally to both SQUIDs, to show the effects of decoherence. The two sets of noise parameters, according 
to Eq[38l correspond to decoherence times of 29 000 and 8 200 time units, giving decoherence times 1/D of 180ns 
and 52ns with our typical time unit. In terms of the formula 1/D = e^R/T, at 20mK these would correspond to 
R « 2Mfl and R « 580fcr2. As expected, the wavefunctions for the states which should not change remained stable 
in these runs. 

In these runs we applied the noise to each SQUID independently. However, one might consider the effects of a 
common noise applied to both SQUIDs together. This would be the case not for the true intrinsic decoherence but 
could represent, say, an instrumental effect due to an external common noise like long range fluctuations in the applied 
field. Interestingly, applying the noise signal so it is the same on both devices seems to have little or no effect. Taking 
the lower curve of Fig. 14 at tsweep = 100 000, one finds that such a common noise produces the same result as the 
independent noise, namely a final overlap of 0.7. This at first surprising result is understandable from the fact that 
the major effect of the noise takes place at level crossing. The noise is mostly effective on the SQUID undergoing 
the sweep and apparently the noise on the control SQUID has little effect. Indeed, turning the noise on SQUID 2 
completely off still leaves the final overlap at 0.7; but turning it off for SQUID 1 leads to a perfect result with overlap 
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FIG. 15: "Good domains" of the {(pT ) plane as in Fig. 11 with the smaller /3 values j3i = (32 = 1.14. Due to the increased 
tunneling the "good" regions where there is a definite state of the SQUIDs have been reduced to small islands. The relative 
narrowness for the diagonal bands, as discussed for Fig. 11, is seen here clearly. For legibility the tableaux are not shown; 
however they correspond to those on Fig. 11. 



SMALLER (3 



As seen in Table II, a reduced (3 gives an increased cotunnei and shorter adiabatic times. Since shorter times gives 
the decoherence less time to act, there is reduced decoherence. This suggest examining the situation with a reduced 
tunnel barrier. We take (3 = 1.14, which according to Table II, increases ojtunnei by a factor six, from 0.0044 to 0.027 
energy units. Fig. 15 shows, analogously to Fig. 11, the good domains of the (^5!^*, (/)2^*)plane for (3 = 1.14. The small 
change in (3 has great effects. The thin black lines of Fig. 11, with "bad wavefunctions ", have now become broad 
bands. Thus starting and finishing points of the sweeps are now restricted to certain "islands". By lowering the tunnel 
barrier, we have created substantial areas represented by the Right panel of Fig. 9, where the wavefunctions arc not 
in a single well with a definite flux in the SQUIDs. While for Fig. 11, this occurred only at the "level crossings" given 
by the black lines, there are now rather large regions with delocalized wavefunctions. Even on the "good" islands 
examination of pictures of the wavefunctions as in Fig. 9 shows that, while they are generally still well concentrated, 
there are little ripples some distance off from the center of localization. Evidently in this parameter range we operate 
on the border to total delocalization. 
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TABLE VI: CNOT sweeps with a shortened adiabatic time, using /3 = 1.14. "Final overlap" is shown for increasing noise 
parameter A. Due to the increased uutunnei permitting faster sweeps, a much stronger decoherence, compared to /3 = 1.19, can 
be tolerated before the sweeps become unsuccessful. Other parameters are as in Fig. 14. The noise frequency was luc ~ 0.042. 
The entry with A = 0.00032 would correspond to a decoherence time of 500 time units, or 3.1 ns with our "typical" parameters. 

To examine the hoped-for improved robustness [3| with respeet to noise/decoherence we show in Table VI the results 

of a series of CNOT sweeps for /3 = 1.14, choosing a sweep on the lower left of Fig. 15 where 

The increased totunnei allows us to set tsweep as short as 1 000 units without violating adiabaticity. Increasing the noise 
parameter A from zero and tabulating the final overlap for state 4, one sees that significantly shorter decoherence 
times as compared to /? = 1.19 (Fig. 14) become possible. The next-to-last entry would correspond to a decoherence 
time of only 500 time units. This reflects in part the phenomenon discussed in connection with Eg 1131 and Table IV 
that the sweep time can be longer than the decoherence time before large decoherence effects occur. 

The 500 time units would be 3.1 ns with our "typical" parameters. This is substantially less than the 20 ns reported 
in ref . [lll | . This is one of our most interesting results concerning the engineering of such devices and perhaps implies 
that the feasibility of such devices is not so very far off. 



CONCLUSIONS 

We have presented a description of one or two interacting rf SQUIDs in a form suitable for numerical simulation 
and a system of programs for carrying out these simulations. In addition to static properties, the behavior under time 
dependent conditions, including flux sweeps for adiabatic gates, and the effects of noise or decoherence, were studied. 

Among the points investigated was the validity of the "spin 1/2 picture" where the behavior governed by the 
full Hamiltonian is approximated by the components of an effective "spin 1/2" system using the lowest levels of the 
double-potential well. The identification between the parameters of the spin Hamiltonian and the full Hamiltonian was 
established and the regime of validity of the approximation were investigated numerically. It is found, using various 
criteria, including "Hilbert space completeness", that the "spin 1/2" picture is approximately valid up to variations 
of the external bias which moves the basic pair of levels a substantial fraction of the principal level splitting. 

For the two SQUID system, we have found it is capable of performing the logical operation CNOT, and have been 
able to determine regions of the ((^f^*, ^2^*)plane with definite configurations of the wavefunctions suitable for CNOT 
sweeps. The conditions for adiabatic behavior under such sweeps were examined and adiabatic times found for typical 
SQUID parameters. The splitting oJtunnei of the basic pair enters into the adiabatic condition and is very sensitive to 
the barrier parameter (3. It is found that the most interesting regions for this parameter are values slightly greater 
than one, and detailed studies were presented for /3 = 1.19. 

Noise or decoherence were simulated as a random flux noise and a formula was derived relating the parameters of this 
random noise to the decoherence parameter D. Numerical simulations verify the validity of the formula. This random 
noise is then applied to the adiabatic logic gates and possible regimes of operation for the gates are found. One of the 
interesting conclusions is that the adiabatic sweep time can be substantially longer than the decoherence time 1/D 
before decoherence effects become very large. This is traced to the fact that decoherence, and also non-adiabaticity, 
are mostly effective at "level crossing" , which is only a small part of the sweep. Indeed for the two-SQUID CNOT 
gate it is found that noise applied to the control bit SQUID, which does not undergo a "level crossing", has almost 
no effect. 

One of our general conclusions therefore, for any devices of this general type, is that "level crossings" should be 
held to a minimum, and when they occur a large tunneling energy or splitting at "crossing" is beneficial both for 
reducing both decoherence and non-adiabaticity effects. 
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It thus appears that an interesting region for the engineering of adiabatic gates involves rather small values of /3, 
where the splitting is large. However the wavefunctions are then nearly delocalized. We have examined the case 
P — 1.14, where the reduced tunnel barrier and thus large tutunnei permit fast operation of the gates. The resulting 
short time for the decoherence to act leads to operations which are successful with presently discussed values of the 
decoherence time, some nanoseconds. However, a careful choice of parameters is necessary. 

Generally speaking, we may say that our studies verify that it is in principal possible to have a regime where one 
can create and manipulate quantum linear combinations of quasi-classical objects, as pictured in Fig. 9. 

While the parameters and examples we have used are specific to the rf SQUID, it will be recognized that our 
methods and results will apply, after some adjustments, to almost all interacting double-potential well systems. 
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